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1 .0 INTRODUCTION 


In order to provide turbulence models useful for computations of the 
flowfields involved in advanced scramjet combustion systems, a number of 
features of these flowfields must be considered. These combustion systems 
involve supersonic flows with embedded subsonic regions and recirculation 
zones, and appropriate turbulence models for scramjet applications must ad- 
dress each of these. The geometry of advanced combustors is often three- 
dimensional, so that the effects of three-dimensionality in the flowfield 
on the turbulence characteristics must be taken into account. Moreover, 
the combustion process in a scramjet system is embedded within a highly tur- 
bulent flow, so that the effects of turbulence on chemical reaction rates 
must be considered; particularly, in the scramjet context, with respect to 
ignition phenomena. On the other hand, to be of maximum utility in scramjet 
combustor design, the turbulence modeling should be as simple and straight- 
forward as is consistent with the requirements of overall accuracy. In this 
application, predictions of mean flowfield structure, the effects of heat 
release, and mean chemical reaction rates are of greatest importance: details 

of the turbulence structure itself can be approximated if the approximations 
introduced do not materially affect the prediction of overall mixing rate, 
chemical reaction rate, and parameters such as the wall skin friction distri- 
bution and flowfield pressure gradient. Since it can be expected that dif- 
ferent effects may dominate in different regions of the flow: nonisotropy 

in recirculation regions; compressibility effects in high speed flow regions; 
and turbulence-chemistry interaction effects in regions in which fuel ignition 
is occurring, a modular approach may be the most efficient turbulence model 
overall. In such an approach, each module contains the turbulence model 
elements which best account for the dominant features of each region of the 
flowfield. 

An assessment of turbulence models for scramjet applications was initi- 
ated in September 1979. During the first year of this work, as outlined 
in Ref. 1, the major effort involved the examination of the multiple dissipa- 
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tion length scale (MDLS) turbulence model, since this approach appeared to 
offer the potential for greater generality than existing models in the context 
of scramjet-related flowfields. In addition to this work, other efforts 
carried out during the first year of this program included the definition 
of a technique for the estimation of the initial conditions required by field- 
equation turbulence models (Ref. 1), an examination of the use of a modified 
dissipation rate equation with the basic k-e two-equation turbulence model, 
and the development of a supersonic-flow compressibility correction to the 
dissipation rate equation in the two-equation (or MDLS) approach. 

Although the results of Ref. 1 indicated that the MDLS model is slightly 
more general than the basic k-e model, the gain did not appear worth the 
added cost of solving two additional equations. Furthermore, the flowfields 
considered in the analyses reported in Ref. 1, while fundamental to and under- 
lying many of the structures found in scramjet flows, did not involve large 
scale recirculation regions. Accordingly, the focus of the second year's 
work shifted to an assessment of the performance of a variety of turbulence 
models in low-speed and high-speed recirculating flows. This work, described 
in Ref. 2, involved the application of several turbulence models to a variety 
of recirculating flows including axisymmetric and planar sudden-expansions. 

The turbulence mdoels studied included the basic two-equation model and the 
algebraic stress model, in both cases with and without modifications to the 
dissipation rate equation designed to enhance the model's sensitivity to 
streamline curvature. Results of this work indicated that the algebraic 
stress model was the most generally applicable in regions of strong recircula- 
tion, and that the modification to the dissipation rate equation, while impor- 
tant in the region surrounding the recirculation zone, had a deleterious 
effect on the overall level of predictions in the region downstream of the 
recirculation zone. On the other hand, the results described in Ref. 2 were 
felt to be inconclusive, since only in planar subsonic recirculating flows 
could clear model differences be discerned. 

A further observation made in the course of carrying out the work des- 
cribed in Ref. 2 was that in complex flowfields it is difficult to separate 
some aspects of the turbulence modeling problem from the numerical problems 
inherent in different computational approaches for solving the governing 
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equations describing the flow. These aspects include the treatment of wall 
boundary conditions, the algorithms used to generate the finite-difference 
form of the equations, and the algorithms used to provide the finite-differ- 
ence solution of the governing equations for the particular turbulence model 
chosen. While not an integral part of the work plan for the program, several 
such problems were encountered and were discussed in Ref. 2. Work in this 
area continued during the program described in this document, with particu- 
lar emphasis on the minimization of numerical diffusion problems in turbulence 
model assessment. 

The work described in this report completes the assessment of turbulence 
models for scramjet applications. Axisymmetric subsonic recirculating flows 
were considered to complete the k-s model and ASM model comparisons initiated 
in the prior program. Swirling flows, while not themselves of direct interest 
in scramjet applications, were found to be stringent tests of turbulence 
modeling so that assessment of the various turbulence models was extended 
to these flowfields. Considerable care was taken to minimize the effects 
of numerical diffusion which, if not carefully minimized, can all but swamp 
differences between turbulence models. Further efforts to incorporate 
higher-order turbulence models in both the time-split and time-unsplit 
MacCormack predictor-corrector schemes were carried out. These efforts 
were unsuccessful : the introduction of stiff equations renders the basic 

solution algorithm unstable, and this result points up the need for turbu- 
lence model development and numerical solution algorithm development to 
proceed in parallel. Just as there is no completely general turbulence 
model, neither does a general numerical solution algorithm exist. Finally, 
methods of modeling scalar transport which do not invoke the Boussinesq 
gradient diffusion hypothesis were investigated. Each of these areas is 
described in detail in this work. 

Although, as noted above, no completely general turbulence model 
exists for scramjet applications, the work described in this report and in 
Refs. 1 and 2 indicates that the algebraic stress formulation is the method 
of choice. This conclusion results from the greater generality of the 
approach, compared with the two-equation model, and from its ability, when 
coupled with appropriate numerical resolution, to model details of the flow 
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such as counter-rotating vortices in a sudden expansion recirculation 
region. In regions of relatively simple flow (i.e. , jets and shear layers) 
the basic two-equation approach remains applicable, and it is easily 
arranged to transition from the ASM to the k-c approach since the turbulence 
kinetic energy and dissipation rate equations are fundamental parts of the 
ASM formulation. Finally, solution of the equation for transport of the 
turbulent species flux, in the same general way as the ASM solution proceeds, 
appears to be a viable method for generalizing these results to more complex 
flows with species and energy transport. However, further model development 
work is requi red in the latter area, and such model development work is 
recommended as an outgrowth of the work described in this report. 
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2.0. ASSESSMENT OF TURBULENCE MODELS FOR SCRAMJET APPLICATIONS 


Three major areas were addressed in this phase of the assessment of 
turbulence models for scramjet applications carried out under this program. 
These areas included assessment of models for subsonic, axi symmetric, swirling 
and nonswirling recirculating flows and the development of the k-e and the 
algebraic stress model approaches for the prediction of supersonic recircu- 
lating flows. The k-e model relates the Reynolds stresses to the mean rate 
of strain through the definition of an isotropic turbulent viscosity, whereas 
the more advanced algebraic stress model calculates the stresses from implicit 
algebraic relationships containing the stresses themselves, the mean rate of 
strain, and the turbulent kinetic energy and its dissipation rate. Modified 
versions of the models employ a new dissipation rate equation whose production 
term was made more sensitive to streamwise curvature effects. A new non- 
equilibrium wall-function treatment was also incorporated into each model. 
These models are discussed in detail in Ref. 2. Swirling flows were investi- 
gated because of the demands they place on models with respect to accurate 
predictions of a wider variety of flowfield features than nonswirling flows 
and their incorporation of multiple shear stress components. In the area of 
supersonic recirculating flows, problems of turbulence model application were 
encountered which highlight the interaction between turbulence modeling and 
numerical solution techniques. 

The assessment of turbulence models for axi symmetric sudden expansion 
flows reported in Ref. 2 was inconclusive since only one diameter ratio was 
calculated using a relatively coarse grid. In this year's program a follow-up 
study covering different diameter ratios and using a finer mesh was conducted 
as the final step in the assessment of turbulence models in subsonic flows. 

Most computational studies to date fail to separate the effects of 
numerics from the influence of turbulence models in assessing their results. 
Current literature therefore abounds with discrepancies in the predictions of 
flowfields with the same turbulence models but different numerical techniques. 
A preliminary study was conducted to address this problem, considering the 
following four areas: 
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1. Selection of test cases 

2. Determination of the solution domain length 

3. Determination of convergence criteria 

4. Reduction of numerical diffusion through grid refinement. 

2.1 SELECTION OF TEST CASES 

Surprising there are very few experimental (or numerical) studies of 
mean and turbulence flowfields in axisymmetric sudden expansions. The only 
applicable study that has detailed velocity and turbulence measurements is 
the Chaturvedi case (Ref. 3). Back and Roschke (Ref. 4) report reattachment 
length measurements in laminar and low Reynolds number turbulent flows 
(^ e ^4000 based on inlet diameter) for 2.6:1 diameter ratio sudden expansions. 
Stuijess and Syed (Ref. 5) provide some centerline velocity data, originally 
obtained by Lipstein and described in a General Electric Co. corporate report 
only. These data are for diameter ratios of 1.6:1, 2.5:1 and 4:1. Finally, 
heat and mass transfer measurements in axisymmetric sudden expansions are 
reported in Refs. 6, 7, and 8. For this work, diameter ratios of 1.33:1, 2:1 
and 3:1 were selected as representative values of the range of diameter ratios 
encountered in practical and research oriented applications. Comparisons with 
experiments are presented for the 2:1 diameter ratio case. 

2.2. DETERMINATION OF THE SOLUTION DOMAIN LENGTH 

Specification of the solution domain length is an important considera- 
tion in computational work. Too short a length causes the exit boundary 
condition to affect the rest of the flowfield more strongly than is physically 
realistic, and an excessively long solution domain results in poor resolution 
or a waste of grid points. Thus a study was conducted to determine a realis- 
tic solution domain length for turbulent axisymmetric sudden expansion flow 
calculations. For this purpose the 2:1 diameter ratio geometry of Chaturvedi 
(Ref. 3) was chosen as the test case, using the k-c turbulence model. This 
grid refinement study consisted of increasing the length of the solution 
domain while maintaining a fixed mesh spacing. More nodes were added in the 
streamwise direction to allow these increases — a necessary step to exclude 
effects due to changes in grid spacing. The value at which the reattachment 
length, the flow parameter chosen to be tested, ceased to vary with further 
increases was taken as the lower limit for the solution domain length. The 
results of this study, presented in Figure 2.1a, indicate a value of 11 step 
heights or more as the recommended length for the k-e model predictions 
of this geometry. A solution domain length of sixteen (16) step heights 
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was adopted for subsequent computations to accommodate anticipated differences 
in reattachment length predictions for the different turbulence models inves- 
tigated. Results of a similar study for the 1.33:1 diameter ratio, given 
in Figure 2.1b, also confirm these observations with respect to solution 
domain length. 

2.3. DETERMINATION OF CONVERGENCE CRITERIA 

In iterative codes the convergence criterion is a measure of the degree 
to which a computed solution satisfies the finite-difference equations. 

For the STEP family of programs (Ref. 9) this criterion is the level of the 
residual sources (SORMAX). A study was conducted, again using the k-e turbu- 
lence model, to determine a realistic convergence criterion for the test cases 
considered in this work. This study was carried out by sequentially decreas- 
ing the level of SORMAX until the change in the reattachment length, the 
flow parameter tested, was less than 1 percent for an order of magnitude 
reduction in SORMAX. For all three area ratios, a satisfactory convergence 
criterion was determined to be 0.001. 


2.4 REDUCTION OF NUMERICAL DIFFUSION THROUGH GRID REFINEMENT 

It is a well -known fact that all upstream or hybrid-upstream* based 
finite differencing techniques, although computationally very stable, intro- 
duce numerical diffusion into the formulation. Unless the effects of this 
artificial viscosity are removed, the solution thus obtained does not satisfy 
the governing partial differential equations and can be in serious error, The 
source of this problem is the truncation error in the one-sided first differ- 
ence used in upstream differencing 
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Note that the leading-order discretization term is equivalent to a physical 


diffusion ter m with an effective diffusion coefficient of 
♦Hybrid. upstream finite. differencing techniques (currently used in many codes 
including the STEP family) use upstream differencing for |Pe ce n| >2 and central 
differencing for -2<Pe C el1i2. Pe ce ll is the cell Peclet number defined in 
Eq. 4. 
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r _ j. uAx , _ 

Numerical ~ 2 

This term is one explanation of the stability of upstream differencing tech- 
niques. By choosing (1) when u > 0 and (2) when u < 0, the discretization 
error terms are always stabilizing. However, they also introduce a numerical 
diffusivity that may dominate the physical diffusion terms. The appropriate 
parameter is the cell Peclet number, defined as the ratio of discretized 
convection terms to discretized diffusion terms 


Pe 


cell 


uAx 

r 


(4) 


where r is the effective (laminar or turbulent) diffusion coefficient. This 
parameter can also be interpreted as 


Pe cel 1 ' 2 


Numerical 


(5) 


Therefore, when Pe ^ = 2, numerical and physical diffusion are of the same 
size. Computationally, if Pe ^ <, 1 the flow is diffusion dominated and 
the form of finite differencing used in the convection (first derivative) 
terms does not affect the results. For |Pe > 2, central differencing 
of the convection terms becomes potentially wiggly, and for convection domi- 
nated flows ( | P e ce i 1 ( > 5), upstream differencing can produce numerical diffu- 
sivities that can interfere with and completely dominate physical diffusion 
terms. 


It is usually quite difficult to accurately quantify numerical viscosity 
effects in multi-dimensional flows. Nevertheless, a good estimate of the 
magnitude of the numerical diffusion coefficient in two-dimensional cases 
can be obtained from Ref. 10 as 


, U R AxAysin20 

Numerical " 777 . 3„ ‘ T77 

4(Aysin 8 + Axcos 8) 


( 6 ) 


where is the resultant velocity, Ax and Ay are the grid spacings in the 
x- and y-directions, respecti vely, and 9 is the angle (between 0° and 90°) 
the velocity vector makes with the x-direction. Several key observations 
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about the nature of numerical diffusion can be made in terms of this expres- 
sion: 

1. Numerical diffusion vanishes when the streamlines are perpendic- 
ular to the mesh, i.e., when the flow is along one of the sets 
of grid lines (6 = 0° or 90°). On the other hand, numerical 
diffusion increases with streamline curvature becoming most seri- 
ous when the velocity vector makes an angle of 45 ° with the grid 

1 ines. 

2. There is no numerical diffusion when the gradient of the dependent 
variable normal to the direction of the flow is zero. 

3. Numerical diffusion gets larger with increasing Reynolds number 
and could become serious in convection dominated flows. 

4. Effects of numerical diffusion can be removed by realigning the 
grid with the flow direction (reducing 9 ), by changing the size 
of the mesh (reducing Ax and Ay) or by using both techniques. 

Recirculating flows are especially susceptible to numerical diffusion even 
if the grid is aligned with the primary flow direction. The presence of 
large regions with strong streamline curvature, such as in the reverse flow 

zone, creates local regions where numerical diffusion effects can be substan- 
ti al . 


Current techniques for the reduction of numerical diffusion involve 
the use of grid refinement (or grid-independency) studies in which the number 
of grid points to be used, to minimize the effects of numerical diffusion, 
is determined separately in each coordinate direction. In this work grid- 
independency was determined in the following manner. First, the number of 
points in the r-direction (NJ) was fixed and the number in the x-direction 
(NI) was increased in increments for a fixed-length computation domain. 

The same procedure was then repeated for a fixed NI and variable NJ. The 
flow parameter tested as a function of grid spacing was the reattachment 
length x R . The number of mesh points for grid-independent results was taken 
as the value of NI and NJ where the reattachment length curve appeared to 
approach an asymptotic limit. A separate grid refinement study of this nature 
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was conducted for the 2:1 and 1.33:1 diameter ratios. Since the 3:1 diameter 
ratio expansion represents a flow in which diffusive effects are smaller 
than in the lower diameter ratio test cases, no separate grid refinement 
study was conducted for this case, and the grid determined for the 2:1 diam- 
eter ratio geometry was used in the definitive calculations of this flow. 

The results of this work are reviewed below. 


2:1 Diameter Ratio Grid Refinement Study 

The results of this study are presented in Figure 2.2. For the x-direc- 
tion, an asymptotic limit for the reattachment length was achieved at values 
of NI equal to 72 or greater as shown in Figure 2.2a. A slightly different 
behavior was observed in the r-direction grid refinement study in Figure - 
2.2b. The reattachment length initially increases with the number of cross- 
stream grid points, reaching a peak at an NJ of 28. The curve then monotoni- 
cally decreases approaching an asymptotic limit at values of NJ exceeding 

52. This behavior is probably due to the non-equilibrium wall-function treat- 

+* 

ment that requires y values of 40 to 100 for plausible near-wall results. 
This range of values for y + is approached only at the higher NJ values as 
shown in Table 2.1. 

TABLE 2.1. Variation of y + Values 


GRID 

SIZE 


NI_ 

NJ 

+ 

y 

22 

22 

250 

22 

24 

218 

22 

28 

181 

22 

32 

161 

22 

36 

144 

22 

42 

125 

22 

52 

103 


Comparison- of velocity profiles provides more insight into the effects 
of numerical diffusion. Figure 2.3 presents the axial velocity profiles 




where k v is the turbulent kinetic energy at the edge of the vis- 


y = - 

cous sublayer, y p is the distance of the node to the wall and v is the kv 
nematic viscosity. 


- 11 - 



Reattachnient length In 9tep heights. 




Figure 2.2. Grid-Independency studies 
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at x/h locations of 2 and 8, as predicted using 22 by 22, 22 by 32, 22 by 
42, and 22 by 52 node grids. The main differences between the profiles are 
in their behavior near the wall, and across the separated shear layer towards 
the centerline. Excluding the coarsest 22 by 22 mesh, the differences between 
the predictions for the remaining grids are small. This is consistent with 
the reattachment length computations reported in Figure 2.2b which also show 
a relatively small variation of x^ with the number of cross-stream grid points 
for meshes finer than 22 by 32. Axial velocity predictions carried out using 
22 by 22, 32 by 22, 42 by 22, and 52 by 22 node grids are reported in Figure 
2.4 at x/h locations of 2 and 8. The differences between these velocity 
profiles are more pronounced than those observed with increasing number of 
grid points in the radial direction since the reattachment length changes 
significantly with the number of streamwise grid points, as shown in Figure 
2.2a. 

This study shows that practically numerical diffusion free and fully- 
converged results can be obtained for a 2:1 diameter ratio axisymmetric sud- 
den expansion, with the STEP code (Ref. 9), by using a 72 by 52 node mesh 
(uniform in each coordinate direction) for a solution domain length 
of 16 step heights and a convergence criterion of a residual source level 
of 0.001. However, it is possible to obtain acceptable results by using 
a smaller number of carefully nonuni formly distributed grid points, as demon- 
strated in Figure 2.1, where a nonuniform 22 by 22 node mesh resulted in 
reattachment length predictions in excess of 8.8 step heights, reasonably 
close to the ultimate value reached in the grid refinement study just des- 
cribed. Thus, the second phase of this grid refinement work involved estab- 
lishing a nonuniform grid which minimized numerical diffusion effects through 
comparison with the results obtained using the 72 by 52 node mesh determined 
previously. Nonuniform 22 by 22 and 59 by 43 node grids were devised by 
studying measurements and coarse mesh calculations of the given geometry to 
identify the regions of sharp gradients, and adjusting grid spacing accord- 
ingly for efficient node distribution. The baseline 72 by 52 node mesh compu- 
tations were then matched with the predictions obtained using these two non- 
uniform grids for both the standard and "modified" versions of the k-e and 
the ASM models. The goals of this study were to document the effects of 
numerical diffusion in terms of reattachment length, velocity, turbulent 
kinetic energy, and Reynolds shear stress predictions. 
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x/h = 8.0 


Expansion, as a Function of 


Calculations of this geometry were carried out for a solution domain 
length of 16.471 step heights, with a residual source convergence criterion 
of 0.001. A sketch of the computational domain is shown in Figure 2.5. 

For all three grids, the same inlet velocity profile extracted from the 
Chaturvedi data (Ref. 3) was used at x/h = -0.471 to start the computations, 
and the inlet turbulent kinetic energy and dissipation rates were calculated 
from this profile using Prandtl's mixing length hypothesis. 

Reattachment Length Predictions 

In Table 2.2 the reattachment length predictions obtained using the 
22 by 22 and 59 by 43 node grids with the 72 by 52 node baseline results for 
all four models are compared. These results indicate that the reattachment 
length predictions with the 22 by 22 grid do indeed suffer from numerical 
diffusion. The reattachment length increases by 16.9%, 16.2%, 14.7%, and 
5.0%, respectively, for the k-e, "modified" k-e, algebraic stress and 
"modified" algebraic stress models when the grid is further refined to the 
baseline 72 by 52 mesh. The relatively small change in the "modified" 
algebraic stress model predictions may be explained in terms of two counter- 
acting phenomena. The fine grid reduces numerical diffusion so that the 
computed separated shear layer spreads less rapidly which in turn produces 
a larger recirculation region. On the other hand, the better resolution 
provided by the fine mesh enhances the sensitivity of the model to the 
secondary strains created by the curvature of the flowfield. Thus, the net 
result is a smaller increase in the reattachment length than that observed 
for the other models. Predictions made with the 59 by 43 node are in excel- 
lent agreement with the baseline results indicating a good node distribution 
within the mesh. The danger associated with not separating the effects of 
numerics from the performance of turbulence models is especially apparent 
for this case. Note that the numerically diffusive 22 by 22 grid predictions 
seem to agree very favorably with the experiment, giving a false impression 
of the predictive capability of the models at this diameter ratio. 

Velocity, Turbulent Kinetic Energy and Reynolds Shear Stress Predictions 

The effects of grid refinement on the velocity, turbulent kinetic energy 
and Reynolds shear stress predictions are discussed below. The results are 
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h = 6.3076 x 10" 2 m 


R 1 = 2h 


Grids : 22 by 22 and 69 by 43 node 
non-uniform meshes, and 72 
by 62 node uniform (In each 
coordinate direction) mesh 

Solution domain length : 16 step heights 

downstream and 0.471 upstream of 
the expansion for a total of 
16.471 stop heights 


R 0 = h 


Convergence criterion: SORMAX i 0.001 


Figure 2.6. Solution domain used In the Chaturvedl (2:1 diameter ratio) predictions 



TABLE 2.3. Reattachment Length Predictions 


MODEL 

GRID 

REATTACHMENT LENGTH IN STEP HEIGHTS 


k-s 

72 x 52 

9.53 



59 x 43 

9.57 



22 x 22 

8.15 


M. k-s 

72 x 52 

9.62 

— 


59 x 43 

9.68 



22 x 22 

8.28 

— 

ASM 

72 x 52 

9.30 



59 x 43 

9.34 



22 x 22 

8.11 

- - 

M. ASM 

72 x 52 

10.22 



59 x 43 

10.21 



22 x 22 

9.73 


Data, Ref. 3 


9.4 + 1* 



* The uncertainty in the Chaturvedi 2:1 diameter ratio reattchment 
length measurements is estimated to be + 1 step height. 
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presented at an x/h location of 5 for the 22 by 22, 59 by 43 and 72 by 52 
node grids. 

Figures 2.6 and 2.7 present the U, k and uv predictions for the k-c 
and "modified" k-e models, respecti vely. Similar results are reported in 
Figures 2.8 and 2.9, respecti vely, for the standard and "modified" versions 
of the algebraic stress model. The following conclusions can be drawn on 
the effects of numerical diffusion in the Chaturvedi predictions from these 
results: 

1. The qualitative effects of numerical diffusion are the same for 
all models. Quantitatively, the "modified" algebraic stress 
model predictions appear to be the least sensitive to grid refine- 
ment. However, this may be the result of the two counteracting 
phenomena discussed in the preceding section. 

2. The agreement between the baseline results and the 59 by 43 grid 
predictions is excellent for all variables and models. 

3. The differences between the 22 by 22 grid results and the baseline 
predictions are significant, and these should be viewed in con- 
junction with the reattachment length computations. 

a. The main effect of numerical diffusion is a more rapid 
predicted rate of spread for the separated shear layer, 
and consequently a shorter reattachment length. 

b. These effects can be seen in the progressively faster rate 
of decay of the centerline velocity. Also higher levels 

of kinetic energy and shear stress are predicted, indicating 
a shorter yet more intense recirculation zone. 

c. A subsequent drop in the levels of these quantities, typical 
for the recovery region, signals reattachment of the separ- 
ated shear layer and beginning of the relaxation regime. 

Summary 

The results of this study showed that the nonuniform 59 by 43 node 
grid does indeed match the baseline predictions extremely well and can be 
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Figure 2.6. U-velocily, K, and uv predictions at X/h = 5.0, k-e 
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Figure 2.7. U-velocily, k, and uv predictions at X/h - 5.0, M. k-c 
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Figure 2.8. U-velocity, k, and uv predictions at X/h = 5.0, ASM. 
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Figure 2.9. U-velocity, k, and uv predictions at X/h = 5.0, M. ASM. 



used in definitive numerical diffusion free predictions of the Chaturvedi 
geometry at substantial savings (%40%) over the 72 by 52 node baseline grid. 
The nonuniform 22 by 22 node grid suffers from the effects of numerical diffu- 
sion, and thus, grids of this coarseness should be used mainly in preliminary 
computations or to establish qualitative flow trends. In addition, this 
study confirmed the role of the reattachment length as a sensitive index 
of the degree of numerical diffusion in both the k-e and algebraic stress 
model predictions. Hence, parametric studies of reattachment lengths as 
a function of grid spacing are sufficient in constructing nonuniform grids 
to minimize numerical diffusion. This important observation was used in 
the subsequent studies. 

1.33:1 Diameter Ratio Grid Refinement Study 

The computational domain length and convergence criteria for this geom- 
etry were established as 16 step heights and SORMAX = 0.001, respecti vely, 
based on the 2:1 diameter ratio results already discussed. The flow geometry 
and solution domain are shown in Figure 2.10. The grid refinement work pro- 
ceeded in the same manner as that already described: the number of points 

in the r-di recti on (NJ) was first fixed, and the number in the x-di recti on 
(NI) was increased in instal lments for the fixed length computational domain 
until further changes in predicted recirculation zone length no longer 
occurred. The same procedure was then repeated for a fixed NI and variable 
NJ. The number of mesh points for grid independent results would have been 
taken as the value of NI and NJ where the reattachment length curve appeared 
to approach an asymptotic limit. However, within the NI (22-78) and NJ 
(22-82) values studied, this limit was not reached. Therefore, as shown in 
Figure 2.11 the reattachment length curves were extrapolated to their 
asymptotic limits using the reported x R values. These results showed that 
a mesh of at least 108 by 130 nodes of uniform spacing was needed to obtain 
grid independent predictions of this geometry. 

The next step was to devise and test nonuniform grids employing fewer, 
but nonuniformly distributed grid points against the reattachment length 
predictions reported in Figure 2.11. Nonuniform meshes of 22, 32, 42, 52, 
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382-0627 8 3-809FR 



Figure 2.10. Solution domain for the 1.33:1 diameter ratio 
axisymmetric sudden-expansion geometry. 



reattachment length in step heights, h X R reattachment length in step heights, 


382 - 062783-809811 



Number of streamwise grid points, Nl 

(A) NJ held fixed at 22 


383-062783-8098R 



Number of cross-stream grid points, NJ 

(B) Nl held fixed at 22 

Figure 2.11. 1.33.1 diameter ratio axisymmetric sudden— expansion streamwise 
(NJ =22) and radial (Nl = 22) direction studies. 
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and 64 nodes in the streamwise direction coupled with 22 uniformly-spaced 
nodes in the radial direction were matched against the asymptotic reattachment 
length predictions in Figure 2.11(a). A similar study, Figure 2.11(b), was 
also conducted using 22, 32, 42, 52, and 52 nonuniformly-spaced nodes in 
the radial direction and 22 nodes with constant spacing in the streamwise 
direction. These tests produced the 64 by 62 nonuniform grid that was used 
in the definitive numerical diffusion free calculations of this geometry 
with the k-e, "modified" k-e, ASM, and "modified" ASM models. 

This work again emphasizes the need for grid refinement and related 
studies before model evaluations are undertaken. The difference in the reat- 
tachment length predictions alone varies by as much as 1.2 step heights be- 
tween the finest (62 by 22) and the coarsest (22 by 22) grids tested. It 
is interesting to note that nonuniform grids simply shift the reattachment 
length curve to higher values at a given value of NI or NJ. The variation 
of reattachment length with grid spacing is remarkably smooth in both coordi- 
nate directions and by and large displays the desired behavior. 

3:1 Diameter Ratio Grid Refinement Study 

Since this diameter ratio represents numerically a less diffusive 
flow than the 2:1 diameter ratio case, the 59 by 43 node grid determined 
for the 2:1 diameter ratio expansion was also used for this diameter ratio. 

The computation domain length and convergence criteria were again taken as 
16 step heights, and 0.001, respecti vely . The flow geometry and the solution 
domain are shown in Figure 2.12. 

Summary 

This phase of the work involved a parametric study to assess the effects 
of numerics on the accuracy of flowfield predictions. Solution domain length, 
convergence criteria and numerical diffusion were' all considered. A study 
of this nature is essential in most computational work to separate the effects 
of numerics from the influence of analytical models in assessing the predic- 
tions. The results of the present study, which provided the solution domain 
length, convergence criteria and the mesh size to be used in the definitive 
test case calculations, are summarized in Table 2.3. 
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Figure 2.12. Solution domain for the 3:1 diameter ratio 
axisymmetric sudden-expansion geometry. 



TABLE 2.3. Summary of Results. 


DIAMETER 

RATIO 

GRID 

SIZE 

SOLUTION 
DOMAIN 
LENGTH, 
Step Height 

CONVERGENCE 

CRITERIA 

1 .33:1 

64 by 62 

16 

0.001 

2:1 

59 by 43 

16 

0.001 

3:1 

59 by 41 

16 

0.001 


It was also shown that the reattachment length is indeed a sensitive index 
of the degree of numerical diffusion for both k-e and algebraic stress model 
predictions. 

2.5 ASSESSMENT OF MODELS FOR AXISYMMETRIC SUDDEN EXPANSION FLOWS AT DIFFER - 
ENT DIAMETER RATIOS 

This phase of the study included the definitive prediction of the 
1.33:1, 2:1, and 3:1 diameter ratio cases with the standard and the "modified" 
versions of the k-e and the algebraic stress models for the computation domain 
length, convergence criteria and the mesh size specified in Section 2.1.1. 
Assessment of model performance and the effects of diameter ratio on sudden 
expansion flow predictions are discussed next in terms of reattachment length, 
mean velocity and Reynolds stress calculations. 

Reattachment Length Predictions 

The reattachment length predictions for the different diameter ratios 
investigated are shown in Table 2.4 and Figure 2.13. Here the "primary" 
recirculation zone is that which forms downstream of the step along the outer 
wall of the sudden expansion; the "secondary" recirculation zone is a counter- 
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TABLE 2.4. Reattachment Length Predictions 


I 

CO 

0 

1 


DIAMETER REATTACHMENT LENGTHS IN STEP HEIGHTS 

RAT! 0 PRIMARY RECIRCULATION ZONE SECONDARY RECIRCULATION ZONE 



KEM 

M.KEM 

ASM 

M.ASM 

MEASUREMENTS 

KEM 

M.KEM 

ASM 

M.ASM 

1.33:1 

7.40 

7.60 

7.73 

8.82 


0.19 

0.18 

0.55 

0.51 

2:1 

9.57 

9.68 

9.34 

10.21 

9.4 + 1 

0 

0 

0.44 

0.41 

3:1 

8.61 

8.66 

8.20 

8.37 


0 

0 

0 

0 



RECIRCULATION ZONE LENGTH/STEP HEIGHT LRZ/h 


15 


O KEM 
□ MKEM 
O ASM 
A MASM 


10 


% 


PRIMARY 

RECIRCULATION 
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DIAMETER RATIO D/d 



SECONDARY 

RECIRCULATION 

REGION 


3.0 


FIGURE 2.13. Turbulence Model Predictions of Primary and 

Secondary Recirculation Zone Length as a Function 
of Diameter Ratio. 


rotating eddy which, with some turbulence models, is predicted to exist in 
the corner between the step and the outer wall. The physical presence of 
this eddy is somewhat controversial ; however certain indirect evidence points 
to its existence. Note from the tabulated data and Figure 2.13 that all 
models predict the existence of a counterrotating eddy at small diameter 
ratios and that the predictions indicate that the size of the eddy decreases 
as the sudden-expansion diameter ratio increases. Prediction of the location 
and size of the counterrotating eddy is a function of both turbulence model 
and grid resolution, providing again evidence of the need for fine grid reso- 
lution in prediction of the details of turbulent reacting flows. 

Note from Figure 2.13 the existence of a maximum in recirculation zone 
length (in step heights) as a function of diameter ratio. This is somewhat 
misleading, since step height also increases with diameter ratio, as h = (D-d)/2. 
In fact, the absolute length of the recirculation zone continues to increase 
as diameter ratio (or step height) increases, as shown in Figure 2.14. A data 
correlation obtained by Drewry (Ref. 11) from a variety of sources is also 
shown on Figure 2.14, and the predictions of all of the turbulence models 
examined are in reasonably good agreement with this correlation. 

Several inferences can be drawn from this comparison of reattachment 
length predictions. First, compared with the differences observed as a 
function of grid spacing with a single turbulence model, the differences in 
zone length prediction between difference turbulence models in this sudden 
expansion flowfield are nearly negligible. Only at the lowest diameter ratio, 
where pressure gradient effects on the mixing process are smallest, do signif- 
icant differences in results of the different models become evident, and 
in this case it is only the modified algebraic stress model which provides 
significantly different results. The difference in results zone length 
between the modified ASM and other turbulence models is consistent with the 
results for planar separated flows described in Ref. 2; however, it must be 
noted that at a diameter ratio of 2.0, where the Chaturvedi data (Ref. 3) 
is available, the recirculation zone length experimentally measured does 
not agree with the modified ASM result (Table 2.4). It might also be 
noted that the computational results (and Chaturvedi's data) are not fully 
in agreement with the correlation band shown in Figure 2.14 and, indeed. 
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FIGURE 2.14. Comparison of Predicted Recirculation Zone Lengths 
with Experimental Data Correlation. 
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the computed results indicate that the linear relationship between recircula- 
tion zone length and step height obtained by Drewry may not, in fact, exist. 

Velocity and Shear Stress Profile Comparisons, Diameter Ratio = 2.0 

Figure 2.15 presents a comparison of the axial velocity profiles as 
predicted by the four turbulence models and Chaturvedi's experimental data 
(Ref. 3) at x/h values of 2, 4, 6, 8, and 12. Up to an x/h of 8, there are 
no significant differences between the profiles as predicted by the four 
models. In fact, referring to the results shown in Section 2.1.1, the differ- 
ences observed throughout the flow are much smaller than differences induced 
by grid refinement. At the final two stations, the standard and "modified" 
versions of the ASM predict larger centerline velocities which correspond 
to a slower rate of spread of the separated shear layer toward the center- 
line. This behavior is more pronounced for the modified ASM, consistent with 
the larger recirculation zone this model predicts. 

Agreement between the measurements and the predictions is generally 
acceptable but by no means perfect in the recirculation and near-center! ine 
regions. Across the separated shear layer, predictions consistently fall short 
of the data indicating a more slowly developing shear layer. KEM and "modified" 
KEM display better agreement with the measurements near and downstream of the 
reattachment point. All models seem to do equally well further upstream. 

Predicted and measured profiles of the axial Reynolds stress component 
are presented in Figure 2.16 for x/h values of 2, 4, 6, 8, and 12. All four 
models compute similar u^u^ profiles along the duct which only differ in mag- 
nitude. Up to an x/h of 8, the KEM, "modified" KEM and ASM predictions are 
hardly distinguishable. In this region the "modified" ASM predictions show 
the most rapid decrease in shear stress level. Beyond x/h = 8, the KEM and 
modified KEM predictions of shear stress decay rate increase and eventually ■ 
these models predict lower shear stress levels than both the modified and 
standard ASM. The behavior of the measured stress profiles is successfully 
predicted by all models. The computed stress levels however are generally 
higher, except at the first station, than the measured values. In the 
recirculation region the "modified" ASM shows the better agreement with the 
data. Further downstream, KEM and "modified" KEM predictions appear to be 
more successful in this respect. 
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FIGURE 2.15. Comparison of Turbulence Model Predictions of Axial Velocity, Diameter 
Ratio 2.0 Sudden Expansion. Data from Ref. 3. 
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FIGURE 2.15, concluded. Comparison of Turbulence Model Predictions of Axial Velocity, 

Diameter Ratio 2.0 Axisymmetric Sudden Expansion. Data from Ref. 3. 
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FIGURE 2.16. Comparison of Turbulence Model Predictions of Axial Reynolds Stress, 
Diameter Ratio 2.0 Symmetric Sudden Expansion. Data from Ref. 3. 
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FIGURE 2.16, concluded. Comparison of Turbulence Model Predictions of Axial Reynolds 

Stress. Diameter Ratio 2.0 Axisymmetric Sudden Expansion. 
Data from Ref. 3. 


Velocity and Shear Stress Profile Comparisons, Diameter Ratio =1.33 


At this diameter ratio, effects of pressure forces can be expected to 
be smaller than at a diameter ratio of 2.0 relative to the details of the 
turbulent mixing process, so that more substantial differences between tur- 
bulence model results than was evident in the diameter ratio 2.0 calculation 
may be expected. 

The axial velocity profiles predicted by the four turbulence models 
considered in this work in the 1.33 diameter ratio case are shown in Figure 
2.17 for x/h values of 2, 4, 6, 8, 10, and 14. There are no experimental 
data available for this diameter ratio. Within the recirculation region, 
excluding the immediate vicinity of the reattachment point, there appears 
to be no significant differences between the profiles as predicted by the 
four models. Beyond the recirculation region, the standard and "modified" 
versions of the ASM predict slightly larger centerline velocities which 
correspond to a slower rate of spread of the separated shear layer toward 
the centerline. This behavior is more pronounced for the "modified" ASM 
consistent with the larger recirculation zone this model predicts. 

Similarly, small differences in model predictions are observed when the 
axial Reynolds shear stress component ITu^ is examined. Figure 2.18. All 
four models predict similar u^uT profiles along the duct which only differ in 
magnitude. In the recirculation region, up to an x/h of 4, the highest stress 
levels are predicted by the KEM models. Here, the "modified" ASM predicts 
approximately 15% lower peak stress levels. Further downstream, the stress 
levels begin to decay. In this region, the standard and "modified" KEM 
models experience a more rapid decay, and eventually predict stress levels 
lower than both the "modified" and standard ASM. These results are also 
consistent with the ITu^ predictions for the 2:1 diameter ratio. 

Results for the diameter ratio 3.0 case are completely consistent with 
those already discussed for diameter ratios of 1.33 and 2.0 and are not 
presented herein. 


- 39 - 



II 




ca * 1 

-.500 0.000 


.500 

U/Uln 


I 

~l 

1.000 


1 

1.500 


2.000 


I 

0.000 


II 

.500 1.000 

U/lMn 


1.500 


2.000 


0.000 


_ 1 r 

.500 1 . 000 

U/Uln 


1.500 


FIGURE 2.17. Comparison of Turbulence Model Predictions of Axial Velocity, Diameter 
Ratio 1.33 Axisymmetric Sudden Expansion. 
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FIGURE 2.17, concluded. Comparison of Turbulence Model Predictions of Axial Velocity, 

Diameter Ratio 1.33 Sudden Expansion. 
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FIGURE 2.18. Comparison of Turbulence Model Predictions of Axial Reynolds Stress 
Component. Diameter Ratio 1.33 Axi symmetric Sudden Expansion. 
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FIGURE 2.18, concluded. Comparison of Turbulence Model Predictions of Axial Reynolds 

Stress Component. Diameter Ratio 1.33 Axi symmetric Sudden 
Expansion. 



Summary 

The results of this assessment can be summarized as follows: 

1. The mean velocity and turbulence field predictions do not change 
significantly with different turbulence models. 

2. In these pressure dominated flows, simple k-e model predic- 
tions appear to be comparable to the more sophisticated ASM 
computations. 

3. At the lower diameter ratio, predictions seem to be more sensi- 
tive to modeling effects; however pressure forces still dominate 
the flowfield. The influence of these pressure forces can be 
appreciated by comparing the axi symmetric and planar reattach- 
ment length predictions at the same effective expansion ratio. 

4’. Performance of the "modified" ASM is consistent in both axi sym- 
metric and planar sudden expansion flows. This model is, how- 
ever, more sensitive to grid refinement in low diameter ratio 
axi symmetric flows. An increase in reattachment length of 15% 
is observed in axi symmetric flows over the standard version of 
the model as compared with a 28% increase in planar backward- 
facing step flows of the same expansion ratio. 
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3.0. ASSESSMENT OF TURBULENCE MODELS FOR SUBSONIC SWIRLING FLOWS 


The results of comparisons of turbulence model predictions in axi sym- 
metric sudden expansion do not produce a clear-cut choice among the turbu- 
lence models tested. This is partially due to the domination of the flow- 
field by the pressure forces, which reduces the sensitivity of the calcula- 
tions to the turbulence models used, and partially to the lack of detailed 
experimental data for model evaluation. A more stringent test of model per- 
formance is the calculation of swirling axi symmetric sudden expansion flows 
where the swirl -induced flow anisotropy combined with the pressure forces 
generated by swirl create a more complicated flowfield for turbulence model 
assessment. While detailed experimental data are again lacking, and swirl 
flows are not of great interest in scramjet applications, assessment of model 
performance in highly swirling flows can provide some information on model 
behavior which can be of use in selecting a generally-useful model for scram- 
jet applications. Thus modifications to the governing equations and the tur- 
bulence models for swirling flows were derived and implemented in the STEP 
codes. Fine grid calculations of a swirling flow of the constant angle type 
for a given swirl number were then carried out using the k- £ and algebraic 
stress models. The selection of the k-e and algebraic stress models for 
this study is significant in that the algebraic stress model represents the 
only general purpose turbulence model outside a full Reynolds stress closure 
(where a differential transport equation has to be solved for each stress 
component) that does not use an isotropic turbulence viscosity concept. 

3 -l • GOVERNING EQUATIONS AND TURBULENCE MODEL FORMULATIONS FOR SWIRLING 

FLOWS 

The governing equations of motion for turbulent axisymmetric swirling 
flows are 

Conservation of Mass 
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where U, V, W are the axial, radial and tangential components of the 
mean velocity, P is the pressure, p and y are the fluid density and dynam- 
ic viscosity, respectively, and u 2 v? w 2 uv, uw, vw are the Reynolds stresses. 
This set of equations, however, is not "closed" due to the appearance of the 
Reynolds stress tensor u..Uj which introduces six additional unknowns to raise 
the total number of variables in the four equations to ten (U.., P, uTuT) . 


- 46 - 


The k-e and the algebraic stress models are used in this context to 
express u..u. in terms of known or calculable variables to "close" these 
equations. 


k-e Model 


The k-e model achieves closure by relating the Reynolds stresses 
to the mean strain rate through the Boussinesq approximation 

3U. 

— | - 2/3 5.J pk 


- P u.u. = 



3x i 


(3) 


The effective (turbulent or eddy) viscosity appearing above, p^, is 
defined in terms of a characteristic length and velocity. If this length 
is taken as the turbulence length scale, k 3/2 /e, and the velocity as k* 5 , p 
can be expressed as 



li t s c^ptfe 

(4) 

where k is the 

turbulence kinetic energy, s is the dissipation 

rate and c 

a constant of proportionality. The Reynolds stresses are then 

defined as 
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for swirling axisymmetric flows. Note that the same single value of u t 
appears in all of the Reynolds stress relationships: this is the isotropic 

viscosity assumption. 
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Algebraic Stress Model 


A complete second-order closure of the Reynolds-averaged governing 
equations would require the solution of a transport equation for 
each of the stress components. Even for two-dimensional flows this can 
be a formidable task, since, in addition to the mean flow equations (equa- 
tions (1) and (2)), eight other transport equations (for IT, 7, w, U7, uw, 
vw, k, and e) need to be solved. Under the assumption that the convection 
and diffusion of each Reynolds stress component can be related to the con- 
vection and diffusion of turbulent kinetic energy, however, the stress 
transport equations can be reduced to a set of implicit algebraic expres- 
sions for the stresses in terms of the mean strain rate, turbulent kinetic 
energy, dissipation rate, and the stresses themselves. The equation for 
the algebraic stress model for a swirling flow, including the near-wall 
correction terms are 


u i u j 
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Equation (11) is the version of the algebraic stress model used in this study. 
The details of the derivation and the underlying assumptions for the model are 
discussed in Ref. 2. 

For axi symmetric swirling flows the algebraic stress model representa- 
tion of the Reynolds stresses become 
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where 


r : s 2/3 Xk/e r = 2/3 pXk 2 /e 

X = (1 - c 2 ) / ( (c 1 - 1) + P/e ) 


a. 

lx 

a 2x 


C 2 C 2 V (1 - c 2> 


C 1 e/k c x / ^ 1 " c 2^ 
k 3 / 2 /e 


K/a 3 / 2 x 


lr 


C 2 C 2 " c 2^ 


i 2r = cj e/k ? r /(l - c 2 ) 


k 3 ' 2 /e 

K/a 3/2 r 


c l> c 2 , Cp and c£ are constants defined in Table 2.6, K is the von Karman 
constant (0.4187) and "a" is the near-wall value of -uv/k (generally taken 
as 0.25). 


TABLE 3.1 Recommended Values for Turbulence Model Constants 


K 

= 0.4187 

c k 

= 0.22 

c 

u 

= 0.09 

C £3 

= 0.36 (c, 

J k 

= 1.00 

C 1 

= 1.8 

J 

£ 

= k 2 /(c - c )c 1/2 

£2 £ 1 u- 

C 2 

= 0.6 

C 

s-1 

= 1.44 

c i 

= 0.5 

C 

^2 

= 1.92 

c 2 

= 0.3 
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k and £ Transport Equations 


Both the k-e eddy viscosity and the algebraic stress model require 
evaluation of the turbulent kinetic energy and its dissipation rate to define 
turbulent time and length scales. The high Reynolds number forms of the k and 
e transport equations used in this study are 

■<* -"*&•>* + 7Jr ) ( 18 ) 

x r 


3e , 3e 

3t + pU 3X 


+ pV 


3e _ 

3r “ p 


c/k (c p - 
£ 1 


£2 


:) + 


3x 




(19) 


where 

k-e Model 


“t + 

\ 3k 
/ 3x 

(20) 

v t . 
^ +u 

\ 3k 
) 3r 

(21) 

a M 

£ 

) Is. 

} 3x 

(22) 

U t + 
£ 

\ 3e 
/ 3r 

(23) 


♦(- £(*))'♦ (&) 2 ] -Wi + H£] (24, 


and 0 ^ are turbulent Prandtl numbers for k and e respectively. These are 
defined in Table 2.6 with model constants c and c . 

ei e 2 
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Algebraic Stress Model 


D k = 
K x 



-t) 

, 3k 

+ ^ 37 

(25) 

D k = 
r 

‘kM 


— 3k 1 
uv f 

ax J 

, 3k 

+ u 3? 

(26) 

D 

£ x 



— 3c 
uv 3r . 

, 3c 

[ + 57 

(27) 

D 

% 

Ul 

[•*£* 

— 3e 
uv 3x 

, 3c 

+ 11 5F 

(28) 




(29) 


and c £ are model constants given in Table 2.6. 


Wall Function Treatment 

Most turbulence models including the present versions of the algebraic 
stress and k-e models are devices for high Reynolds number flows. However, 
in the vicinity of solid boundaries where the velocities are small, the low 
Reynolds number effects previously neglected become significant and should be 
accounted for. This can be accomplished either by solving the low Reynolds 
number form of the transport equations or by developing wall functions that 
introduce these effects into the existing high Reynolds number models. 

Chieng and Launder (Ref. 12) found that the first option required vast amounts 
of computer time due to the slow convergence characteristics of the low 
Reynolds number models. On the other hand a new wall function treatment pro- 
posed by the same authors was shown to incorporate these effects with prac- 
tically no increase in computing time. An expanded version of this treatment 
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is used in the present study. Details of this approach are given in Ref. 3. 

3.2 MODEL ASSESSMENT IN SUBSONIC SWIRLING FLOWS 

Model assessment for swirling flow computations was carried out through 
fine grid calculations of a constant angle swirl flow at a swirl number of 2. 


Swirl ing 

flows 

are usually characterized by defining a swirl number 

which is 


the ratio 

of the axial flux of tangential momentum to the axial flux 

of axial 



momentum 

times 

the inlet radius. The definition of the swirl number 

is then 




6. 





c - <t> 

- G x R i 

(30) 


where 




— 

G d> 

= 

f R 1 (Wr) pU 2:rr dr 
•'R 


— 

G x 

= 

f R i f R i 

/ (Up) U 2irr dr +/ 1 p 2irr dr 

•'r j R 


— 

U 

= 

Axial velocity component 



W 

= 

Tangential velocity component 



R 

= 

Inner radius (hubbed swirlers) 



R. 

i 

= 

Inlet radius 



P 

= 

Local density 



P 

= 

Local static pressure 


— 


A number of flow phenomena are related to and affected by the presence of 
strong swirling in combustor flows. These phenomena are complexly inter- 
related and can have profound effects on performance parameters such as 
combustion efficiency, pressure losses, flammability limits, combustion flow 
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stability, and nozzle thrust losses. While a swirl number of 2 is extremely 
high from the standpoint of practical combustor flows, assessment of turbu- 
lence models in high swirl number non-reacting swirling flows provides a use- 
ful step in developing turbulence modeling for complex flows, for at high 
swirl numbers the differences between models are maximized. This suggests 
the utility of experiments run at swirl numbers that would otherwise be of 
little interest from a propulsion system standpoint. 

Two-equation (k-c) model and algebraic stress model calculations were 
carried out for a uniform density flow in a 2:1 diameter ratio sudden expan- 
sion with and without swirl. The computational domain and flow parameters 
for these calculations were as shown in Figure 2.19: inlet swirl was of the 

constant-angle type. No separate grid-refinement tests were conducted for 
these calculations; however, based on the criteria given in Section 2.1.1, 
these results should not be significantly affected by numerical diffusion. 

The effects of swirl on the overall flowfield at this swirl number 
is demonstrated by the streamfunction plots given in Figures 2.20 and 2.21. 

The formation of a large centerline recirculation zone and the enhanced radial 
mixing shown are both direct consequences of the axial momentum deficit crea- 
ted by the non-zero tangential component of the mean velocity. Significant 
differences between the k-e and the algebraic stress models appear only in 
the swirling flow calculations. This supports the hypothesis noted earlier 
that the two models would show a different sensitivity to flow anisotropy, 
induced, in this case, by swirl. While detailed experimental data does not 
exist for this configuration, experiments such as those discussed in Ref. 13 
indicate very large centerline recirculation regions at high swirl numbers, 
indirectly supporting the ASM prediction of the swirling flow. It might also 
be noted that a large recirculation region along the centerline such as seen 
in Figure 2.21 would be unstable to small disturbances (and thus would be 
subject to axial position fluctuations) (Ref. 13). Note further that near 
the wall the predicted flowfields are very similar, so that measurements of, 
for example, static pressure distribution would not necessarily reflect the 
large centerline differences between these flows. Thus detailed measurements 
within the flowfield will be required to provide data with which to directly 
test the models. 
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R = h 


GRID Nonuniform 61x42 node 

Swirl Number = 2.0 

Swirl Severity Parameter* 

= 1.0 

Reynolds Number, Re. 

= 1.0x105 h 

Inlet Turbulent Intensity 
= 0 . 1 % 


FIGURE 3.1. Computation Domain and Flow Parameter. 


Swirl severity parameter is the ratio of the maximum tangential velocity 
component to the maximum axial velocity component. 
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The axial velocity profiles shown by Figure 3.4 compare the swirl and 
the no swirl cases for the two models. The k-e model predictions indicate a 
large (4.68 step heights long) irregularly shaped centerline recirculation 
region shown in Figue 3.2 that extends well into the inlet and appears to 
signal an eventual breakdown on the existing vortex at higher swirl inten- 
sities. The algebraic stress model calculations on the other hand show an 
actual breakdown of the centerline recirculation region into three discrete 
vortices — two small eddies appear in the inlet and a large primary eddy 
(6.37 step heights long) is positioned roughly 0.76 step height downstream 
of the expansion plane (Figure 3.3). The tangential velocity profiles pre- 
sented in Figure 3.5 show the same qualitative behavior for both models, 
however, the magnitude of the velocities are significantly lower for the 
algebraic stress model. The turbulent kinetic enery distributions plotted 
in Figure 3.6 also show roughly the same behavior for both models; here, 
however, the algebraic stress model predicts significantly higher levels of 
turbulence intensity along the combustor axis. The Reynolds shear stress 
profiles given in Figure 3.7 appear quite similar with the k-e model, pre- 
dicting slightly higher stress levels in the near-field flow and the algebraic 
stress model further downstream. These predictions are consistent with the 
flow pattern depicted in Figures 3.2 and 3.3. 

Summary 

The computations just described show that significant differences be- 
tween the k-e and ASM models exist in flows in which anisotropic viscosity 
effects can be expected. For flows in which there is only one major shear 
stress component, such as the conventional sudden expansion in regions away 
from the recirculation zone the additional complexity involved in the ASM 
formulation may not be warranted. Note, however, that even in simpler flows 
there are regions in which secondary shear stress components may be nonnegli- 
gible. Just such a region is the near-corner region, and it will be recalled 
that in this region the ASM formulation indicates a larger secondary eddy 
than does the k-e approach. 
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(b) ASM Predictions 


FIGURE 3.4. U-Velocity Profiles 
















4.0. ASSESSMENT OF TURBULENCE MODELS FOR SUPERSONIC 
RECIRCULATING FLOWS 


Previous efforts (Ref. 2) to incorporate the k-e model in to the 
TWODLE code of J.P. Drummond (Ref. 14) for supersonic recirculating flow 
calculations were hindered by numerical stability problems. Eventually, 
a scheme was devised that appeared to be stable but required very small 
time steps due to the explicit nature of the technique. In addition, the 
fact that the k and e transport equations were solved in the physical co- 
ordinates introduced inaccuracies for irregular geometries. Two options 
were available: modify this scheme so that the k and e equations are 

solved implicity and simultaneously in the transformed coordinates, or 
implement the k-e model in a new implicit version of the TWODLE code (Ref. 
14) that does not use time-split finite differencing. The second option 
was investigated as part of this year's work. 

4.1. GOVERNING EQUATIONS AND WALL FUNCTIONS 

Mean Flow and Turbulence Equations 

The mean flow and turbulence model equations for the k-e model can be 
written as follows for two-dimensional elliptic planar flows: 

Conservation of Mass 

It + i7 oU + £ (3D 


x-Momentum 


Ft oU + 3x ° m 

- - H + FT [ 2 “T f - 2/3 »t (f + I) - ^ 
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y-Momentum 


3pV , 3 mi, . 3 i,2 

IF 3x pUV + 57 pV 


3P + _l|\. + iv\l j_ 3V 

3y 3x [ y T \3y 3x jj 3y [ 2y T 3y " 


(33) 

2/3|j T (lx' + ly ) ' 2/3pk ] 


k Transport Equation 


+ pVk 


pP - P£ + 




(34) 


e Transport Equation 


3ps , 3 I, , 3 ,, 

TT + 3x pU£ + 37 pVs 


C p 

C p^-P-C + 

£i k e 2 k. 


w'[(!r +p )§] + w[(57 + 1I )f7] 


(35) 


where 

U 

V 

P 

k 

£ 

P 

y 

p 


streamwise mean velocity component 
transverse mean velocity 
pressure 

turbulence kinetic energy 

turbulence kinetic energy dissipation rate 

density 

dynamic viscosity r - ? jn 5 

production rate of k, ^ {z[(f ) + ($/]+(§ + f ) 

-! (i + |) 2 | -^(IM) 
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Uj = u t + u (total viscosity) 
k 2 

y t = co — (turbulent viscosity) 

and a £ are the Prandtl numbers for k and e, respectively, and c , 
c and c, , are constants. 1 

Following Ref. 12 these equations can be put in TWODLE form by defining the 
ft, F, (a and ft vectors as 


31) 3F 3(1 _ 

3t 3x 3y " 

where 
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pu 
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pv 

PUV + Ty X 
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y 
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P (p - e) 

pe/k(c P - c e) 


(36) 
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P + 2/3u t D - 2 u t + 2/3pk 

r = . /1U + 1V\ 

T yx " U T \ 3y 3x / 


P + 2/3y T D - 



U £ 


2y T + 2/3pk 
1 ay 
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3x 


3e 
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3k 
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3e 
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c 


u 


c ki 


(c„ - c c )c 

c2 £ 1 ; 


1/2 


311 . 3V .. . 

gj + compressible flows 

0 incompressible flows 


D 2 - 2/3kD 


Currently recommended values for the constants are 


c = 

u 

0.09 

\ = 

0.22 

c = 

£i 

1.44 

O 

fT) 

II 

1.92 

K 

0.4187 


with these values <jl, and <j become 

K £ 

a k f = 1.0 

a * 1.217 

£ 


Wall Functions 


The near-wall low Reynolds number effects are incorporated to the k-e 
model using the non-equilibrium wall functions of Chieng and Launder (Ref. 12). 
However, since TWODLE uses nodal values rather than control volume averages ,^ 
some changes were made in implementing these wall functions. A typical near 
wall region is shown in Figure 4.1. Here it is assumed that node w is at 
the wall, and w + 1 is in the fully turbulent region 


^w+1 



v 


20 


(37) 


where y w+1 is the distance from node w + 1 to the wall and k y is the turbulent 
kinetic energy at the edge of the viscous-sublayer. The level of k y is obtained 

by extrapolating the line through k w+1 and k w+2 to y=y v , hence 
~z 

1 actually becomes 0.614, however a value of 1 is more common in the 

K literature. There are no significant differences in the predictions 
obtained with these two values. 

Chieng and Launder wall functions use near-wall cell integration to cal- 
culate mean production and dissipation rates for those cells. 
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(38) 


k 


v 


k + y w+l " y v 

W+1 V2- Vl 


(k, 


w+2 


" k w+l^ 


The thickness of the viscous sublayer, y v> and the mean flow velocity at that 


location, U v , were then calculated from 


y v = v Re v /k v » 


(39) 


and 


U v = Re v (V p)/ C 2 * 


(40) 


Re/ is assumed to be 20, and the wall shear stress x w is given by 

T w = K * 0 Vl k v H/ ( U E * Vl kv 


(41) 


where K = 0.4187 c^ 2 and E = e k Rev /Re v . The production rate of k at w+1 
is then calculated as 


w+1 




U - U 
w+1 v 


(42) 


W+1 \Vl " y v, 

and the near-wall dissipation rates are expressed, following Spalding (Ref. 15) 
and Pope and Whitelaw (Ref. 16), respectively, as 


'w+1 


^ 3 / 2 /=.y, 


'w+1 


Sr W+1 


(43) 


and 


'w 


2 v 


(44) 


T y k ^ 

The universal viscous sublayer thickness constant, v 
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4.2 MODEL IMPLEMENTATION AND TEST CALCULATIONS 


Model Implementation 

The k-s model and the associated wall functions discussed in Section 
2.3.1 were implemented into the solution procedure of the code in such a 
way as to set up a sequential calculation of the velocity and turbulence 
fields. In this scheme, the density and velocity fields are obtained 
first. The turbulent kinetic energy equation is solved next in the same 
way as for the velocity field. Then the dissipation rate equation is in- 
tegrated, using the latest density, velocity and turbulent kinetic energy 
fields. Finally, the turbulent viscosity is updated and the source and 
sink terms for the k and e equations are calculated. This scheme intro- 
duces three additional subroutines to the code: 


1. MUKEM: Calculates the turbulent viscosity from its 

2 

definition, y = c p k /z 

U 1 1 

2. PRODK: Calculates the production rate of turbulent 

kinetic energy, 

? = - U t /P 

3. SOURCE: Calculates the source-sink terms for the k 

and e transport equations, p(P-e) and 
p s/k (c p - c^e), respectively . 



The solution procedure thus becomes 


CALL INTEG (ISS=1) 

CALL INTEG (ISS=2) 

CALL INTEG (ISS=3) 

CALL MUKEM 
CALL PRODK 
CALL SOURCE 

Advance the time and repeat. 


Solve for velocity and density fields 

Solve for k field 

Solve for e field 

Calculate u t 

Calculate P 

Calculate p(P-e) and p e/k (c P - c e) 

£i £2 
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There are, of course, many variations within any given scheme that 
can be tried to improve stability and/or numerical efficiency. This parti- 
cular scheme was chosen for its proven stability characteristics in ellip- 
tic flow calculations with other computational techniques. The flow chart 
of the TWODLE code with the k-e model is given in Figure 4.2. 


Test Calculations 

Initial testing entailed assessing the compatibility and the logic 
of the changes to the code with and without the k-e model. The first step 
was to reproduce, with one of the existing algebraic viscosity models, pre- 
vious predictions of the Mach 5 10° compression corner test case. Upon 
completion of this task, detailed assessments of the k-e model in this 
application were to be undertaken. Preliminary calculations of the Mach 5 10° 
compression corner test case with this model showed severe stability prob- 
lems. To track down the cause of the problem, it was decided to simplify 
the flow geometry further, to a Mach 5 flow between two parallel plates, 
and specifying a uniform grid by bypassing the coordinate transformation. 

This geometry represents possibly the simplest internal flow test case 
that still contains the features needed to evaluate the k-e model. The 
use of a uniform grid eliminates any potential problems due to coordinate 
transformation and also simplifies the computation of the source-sink terms. 
Calculations of this geometry were done in stages to isolate and identify 
the source of the stability problems. The k and e transport equations can 
be written in the following general form 


9p4> . 

at 

i 


PU |i ♦ 




II 


V 


III 



(45) 


where $ 


D 




k or e 
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FIGURE 4.2 TWODLE Flow Chart. 
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Numerically, the transient, I, convective, II, and diffusive transport, 

III, terms are stabilizing in nature, whereas the source terms, IV, 
especially when stiff, can cause stability problems. The first step in 
the investigation was to set the source terms identically to zero and 
carry out the calculations with only the transient, and convective and 
diffusive transport terms activated. These calculations produced stable 
results for this initial/boundary value problem. The next step involved 
predictions with the complete transport equations for k and e except that 
the production term, p, in the source terms was set to zero. These re- 
sults also proved to be stable. The final test was calculations with the 
complete k and e transport equations. These calculations, however, experi- 
enced stability problems almost immediately, leading to negative k and z 
values and to eventual collapse of the solution scheme within 120 time 
steps. This clearly shows the adverse effects of the stiff source terms 
on the solution scheme. Changes in the wall boundary treatments, and vari- 
ations in the relative order of calculations of the source terms within the 
solution scheme failed to improve the outcome. 
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These results, combined with previous experience, indicate severe 
stability problems in the solution of source-dominated (stiff) transport 
equations with both time-split and time-unsplit MacCormack predictor- 
corrector algorithms. This is a serious problem since most advanced tur- 
bulence (k-e and higher-order closure schemes) and combustion kinetics 
models inherently contain stiff partial differential equations. The ef- 
forts in this respect described in this report were directed to investi- 
gate and devise stable implementation techniques for the k-s model and its 
boundary conditions within the given solution scheme of the code. In other 
words, the goal was to tailor the model to fit the code. These efforts 
have failed due the inability of the solution algorithm to negotiate, in 
any direct way, the effects of the stiff source terms in the k and c (or 
in any other source dominated) transport equations. Other researchers 
have ail so encountered similar stability problems with MacCormack 
predictor-corrector type block solution algorithms. These findings sup- 
port the results described in Ref. 2, in which a compromise scheme that 
solves for the mean flow and the density fields with the MacCormack tech- 
nique, and then evaluates the k and e equations explicitly was recommended. 
Overall, this study raises two fundamental questions. One is the apparent 
lack of versatility of this widely used computational technique which seems 
to be severly restricted in its applicability to advanced turbulent fluid 
flow analysis. The second question is the relative merits of the current 
practice of separating the development of physical models from research 
in numerical analysis as opposed to a unified approach that closely coordi- 
nates developments in computational techniques with the advances in physi- 
cal models to ensure compatibility in their eventual application to scien- 
tific and engineering problems. 
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5.0. ASSESSMENT OF MODELS FOR TURBULENT SCALAR TRANSPORT 


Transport of a scalar mean quantity C (heat, species, etc) is described 
by the equation 


3oC 

3t 


sU i 


3C 

3x i 


3 

3x. 

i 



°v ) 


(46) 


where U^. is the mean flow velocity, p is the density, Y is a molecular trans- 
port coefficient, and u^c is the scalar flux correlation which is the turbu- 
lent scalar transport counterpart of the uTu” correlation of the momentum 
equations. Closure of the above equation requires either solution of a 
transport equation for u^c or the modeling of this quantity in terms of 
known or calculable variables. 


5.1. CLOSURE OF THE SCALAR TRANSPORT EQUATION 

A transport equation for the scalar flux uTc can be derived by multi- 
plying the equation for the instantaneous value of the scalar C (=C+c) by 
U^ and adding it to the x. -component of the Navier-Stokes equations multi- 
plied by c. Upon ensemble averaging, the results may be expressed as 
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( 47 ) 


- sf- ( u.u.c + 6 . . 

oXj \ i J P 


ij ) 


VI. 

Diffusion Transport 


This equation is valid for incompressible flows and where gradients in C 
are small enough for Y'/Y and v‘/v to be entirely unimportant (Y' and v' 
being the fluctuating parts of Y and v) and for p'/p to be significant 
only in the gravitational term. The first and the second terms in this 
equation are exact and require no modeling. The buoyant generation term, 
III, is conveniently modified as follows 

(48) 

where the dimensionless coefficient a is defined as* 

C 


The dissipation term, IV, is zero in isotropic turbulence and is negli- 
gible in non-isotropic turbulence provided that the turbulence Reynolds 
number is high. The pressure-scalar-gradient correlation, V, is the 
counterpart of the pressure-strain correlation in the stress equations. 

With direct dissipation negligible, this provides the mechanism which limits 
the growth of fluxes. Finally, term VI denotes the rate of spatial trans- 
port of u.c due to velocity and pressure fluctuations. This level of clo- 
sure corresponds to a full Reynolds stress formulation for momentum trans- 
port and requires modeling of the pressure-scalar-gradient correlation, V, 
and spatial transport, IV, terms. Modeling of these terms are discussed in 
some detail in Refs. 17 and 18. 


3p 

"C 



★ 

If C stands for temperature and the fluid is an ideal gas, a is unity. 
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Solution of transport equations for scalar flux components is both 
conceptually possible and economically feasible for most flowfields. Simi- 
lar levels of closure however also exist and are widely used in computa- 
tional work. These are discussed next. 


5-2. ALGEBRAIC STRESS MODELS FOR SCALAR TRANSPORT 

These schemes, the counterpart of the hydrodynamic algebraic stress 
models, employ algebraic stress modeling rather than a complete differen- 
tial closure of equation (47). These models for the scalar fluxes, given 
here for the temperature field, generally have the following form (Ref. 18) 
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This form of the model requires solution of transport equations for 
T ' 2 and £y which can be derived in a similar fashion to the k and s equa- 
tions (Ref. 1 ). Further simplifications however are possible by relating 
£y to e, k and T ' 2 by 

£ T = cT - £/kT ' 2 (50) 
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and approximating T' 2 , when the dissipation and production rates of T' 2 
are nearly in balance, as 


T-2 




3T 

3x k 


(51) 


The recommended values for the model constants Cj , c 2 and Cj' are 3.2, 
0.50, and 1.6, respectively. This model, first p?opos5d 7 years ago, has 
yet to be tested in recirculating flows. Free shear flow predictions with 
this model (Ref. 18) however, exhibit the correct behavior of progressive 
collapse for horizontal buoyant mixing layers and surface jets as the mean 
Richardson number increases, and also show reasonably good agreement with 
published turbulent wake and plane jet thermal data. These models are yet 
to be applied to recirculating flow calculations. Further model develop- 
ment, refinement and testing is needed to explore the full potential of 
this level of closure. 


5.3. GRADIENT DIFFUSION MODELS 

These simpler models relate the scalar fluxes to the hydrodynamic 
turbulence properties and the mean scalar gradient 


u .c = u .u . — 
i i j e 3x. 

J 


(52) 


or in terms of an effective turbulent Prandtl number, G 
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- u i c 
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2 

where v = c k /e 
t u 


This is the scalar- transport model that nearly all 
practical calculation schemes have adopted to date. Experiments suggest 
G.j. is approximately two-thirds in many free shear flows but is some 50% 
higher in the vicinity of a wall. The biggest drawback of this model is 


that it is based, like the k-e model, on an isotropic turbulent transport 
coefficient concept and thus is theoretically limited to free shear flow 
and some boundary layer calculations. 
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6':0 SUMMARY: TURBULENCE MODELS FOR SCRAMJET FLOWFIELDS 


The overall objective of the work described in Refs. 1 and 2 and in 
this report has been to establish appropriate turbulence models to use in 
computations of scramjet flowfields. These combustors involve three-dimen- 
sional reacting flows with embedded recirculation regions, and the fuel-air 
mixing rate is critical to the overall performance of the combustor. Thus 
to be usable in design evaluation and data interpretation, turbulence models 
must be reasonably accurate over a broad range of conditions. 

Based on the work reported in Refs. 1 and 2 and in this report, the 
recommended turbulence modeling approach for use in scramjet calculations 
is the algebraic Reynolds stress model. In regions of strong streamline 
curvature, this model should be used with the dissipation equation modifica- 
tion described in Ref. 2 which were designed to improve the sensitivity of 
the model to these effects. The algebraic Reynolds stress approach is par- 
ticularly valuable where multiple stress components are important, such as 
in three-dimensional flows, since the basic two-equation k-e model involves 
an assumption of effective viscosity isotropy that is not borne out by ex- 
perimental results. Where only a single stress component is non-negligible, 
the two-equation approach provides good results. Since the algebraic Rey- 
nolds stress formulation involves the solution of the turbulent kinetic 
energy and dissipation rate equations, it is easily arranged to allow the 
algebraic stress formulations to relax to a two-equation model as the dif- 
ferent stress components become negligible in a given calculation. 

This work has also indicated the extreme importance of considering the 
interaction of the turbulence model and the numerical solution procedure as 
key parts of the development of scramjet combustor models. The interaction 
takes two forms: compatibility of the turbulence model and the solution 
procedure, and the interaction of the model -predicted diffusion with that 
generated artificially by the solution procedure itself. Compatibility 
issues arise because the models recommended in this study involves the solu- 
tion of source-dominated, "stiff" transport equations. The numerical solu- 
tion procedure must be able to accept stiff equations to be compatible with 


- 80 - 



these models. Numerical diffusion is a feature of most numerical solution 
procedures and under some circumstances it can dominate the diffusion pre- 
dicted by the turbulence model. This must be avoided by careful attention 
to grid size and location. Finally, establishment of the proper initial 
and wall boundary conditions is critical to proper use of turbulence models. 
This has been emphasized throughout this work, and the techniques reported 
herein and in Refs. 1 and 2 for establishing initial and boundary conditions 
are recommended for future scramjet modeling work. 
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7.0. CONCLUSIONS AND RECOMMENDATIONS 


The conclusions reached as a result of the work outlined in this 
report can be summarized as follows: 

7.1. SUBSONIC AXISYMMETRIC RECIRCULATING FLOWS 

1. For diameter ratios of 2:1 and 3:1, in axisymmetric flows, 
mean velocity and turbulence field predictions do not 
change significantly with different turbulence models. 

This can be ascribed to the dominance of pressure forces 
over turbulence diffusion in these flowfields. 

2. Axisymmetric sudden expansion flowfield predictions are 
especially sensitive to grid refinement. The numerical 
diffusion inherent in a coarse grid can produce effects 
on the predicted overall mixing process which are larger 
than those produced by variations in turbulence models. 

3. The algebraic stress model produces results in axisymme- 
tric flows at low diameter ratios (where pressure effects 
are reduced) which are consistent with its performance in 
planar backward-facing step flows. Thus the ASM with modi- 
fications introduced to increase the model's sensitivity 

to streamline curvature can be recommended in regions 
where strong recirculations exist, and the unmodified ASM 
in other regions. 

7.2. SUBSONIC SWIRLING FLOWS 

1. Although of little direct interest in scramjet applications, 
swirling flows at large swirl numbers are stringent tests 
of modeling: large differences between k-e model and ASM 

model predictions are evident. This is due to the strong 
anisotropy of the stress components in swirling flows, which 


- 82 - 



is not adequately accounted for by the effective viscosity 
assumptions inherent in the k-e model. 

2. Comparison of model predictions with experimental results 
for strongly swirling flow could provide the means to fur- 
ther develop and improve the ASM model. Appropriate exper- 
imental data is not now available. The ability of the ASM 
to handle strongly nonisotropic flowfields is of potential 
value in the modeling of 3D flows in scramjet combustors, 
even in the absence of swirl. 

7.3 SUPERSONIC RECIRCULATING FLOWS 

1. Efforts, described in this report, to devise stable tech- 
niques for the implementation of the k-e model and its 
boundary conditions in both the time-split and time-unsplit 
MacCormack predictor-corrector algorithms were unsuccessful. 
This result is apparently caused by the inability of the 
solution algorithm to negotiate, in any direct way, the 
effects of stiff source terms in the k and e transport 
equations. 

2. A more detailed study is needed to identify the required 
changes in the MacCormack solution scheme that will elimi- 
nate the drawbacks identified in this work. The present 
work raises fundamental questions with respect to the merits 
of separating the development of physical flowfield models 
from research in numerical analysis. Coordination is re- 
quired to ensure that numerical schemes and physical models 
are ultimately compatible. 

7.4. MODELING OF TURBULENT SCALAR TRANSPORT 

1. Solution of the transport equations for scalar flux compo- 
nents is both conceptually possible and economically feasi- 
ble. Schemes that employ the algebraic stress closure 
rather than a complete differential closure are appealing 
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in terms of consistency with the level of closure chosen 
for the hydrodynamic field. 

2. Further work is required to assess the potential of alge- 
braic scalar flux modeling in more complicated flov/fields 
such as recirculating sudden-expansion flows. 
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